#!/usr/bin/env Rscript

# This is R code for project nonccRCC_391
# Load common functions and libraries
source("/gnet/is3/research/data/dnaseq/analysis/thongn/projects/salthill/Rlib.R")
cmd_args = commandArgs(TRUE)


## ======================================================================
## DN/DS calculation
## ======================================================================
# The method for this calculation is from Ostrow et al PLoS Genetics 2014.
# Briefly, I calculated number of synonymous (s) and nonsynonymous (n) mutations, and number of synonymous sites (S) and nonsynonymous sites (N) across all genes containing somatic mutations.
# Then dN/dS is calculated as (n/N)/(s/S).

load("/gnet/is3/research/data/dnaseq/analysis/thongn/projects/nccRCC_391/dnds/dnds.rda")
barplot(omega, names.arg = subtypes, border = NA, col = cols, ylim = c(0, max(omega)*1.05),ylab = "dN/dS")


## ======================================================================
## Mutation rate
## ======================================================================
load("/gnet/is3/research/data/dnaseq/analysis/thongn/projects/nccRCC_391/rate/muta_rate.rda")

qp <- ggplot(rate_final, aes(x=Subtype2, y=rate, fill=Subtype2)) + 
    geom_boxplot(outlier.size = 1.5) + 
    # geom_jitter(size = 1.5) + 
    xlab("") + ylab("Mutation rate (mutations per Mb)") + theme_bw() + 
    theme(axis.text.x=element_text(angle=60, hjust=0.5, vjust=0.5)) +
    guides(fill=FALSE) +
    facet_grid(. ~ type) +
    scale_y_continuous(limits=c(0,11), breaks=seq(0, 11, 1))
print(qp)



## ======================================================================
## Simulation
## ======================================================================
load("/gnet/is3/research/data/dnaseq/analysis/thongn/projects/nccRCC_391/simulation/simulate.final.rda")
names(simu_list) = c("ccRCC (TCGA)",  "pRCC", "chRCC", "RO" )
cols = c("orange",  "violet", "red", "green")
get_lim <-function() {
    x = c(unlist(simu_list), observed_freq)
    c(min(x), max(x))
}
boxplot(simu_list, vertical = TRUE, outline = FALSE, border = cols , ylab = "", ylim = get_lim())
stripchart(simu_list, method = "jitter", pch = 20, cex = 0.4, col = cols, vertical = TRUE, add = TRUE, jitter = 0.2)
points(x = 1:4, y= observed_freq, col = "black", pch = "*", cex = 2.5)
mtext("Frequency of Deleterious Mutations", side = 2, line = 2.8)



## ======================================================================
## Coverage plot
## ======================================================================

# Load sample information from sup_table to get subtype information
sam_info = read.xls("/gnet/is3/research/data/dnaseq/analysis/thongn/projects/nccRCC_391/04.nccRCC-SupTables.08122014-v3.xls", sheet = 3, skip = 1)
sam_info2 = sam_info[, c("Exome.Tumor.Sample.ID","Subtype")]
colnames(sam_info2)[1] = c("SamId")

coverage = read.xls("/gnet/is3/research/data/dnaseq/analysis/thongn/projects/nccRCC_391/04.nccRCC-SupTables.08122014-v3.xls", sheet = 4, skip = 1) 
coverage = coverage[, c("Sample.ID", "Mean.coverage")]
colnames(coverage) = c("SamId", "cov")

cov_df = merge(coverage, sam_info2, by = "SamId")
cov_df = data.frame(cov_df, all = "All samples")

qp <- ggplot(cov_df, aes(x=Subtype, y=cov, fill=Subtype)) + 
    geom_boxplot(outlier.size = 0) + 
    geom_jitter(size = 1.5) + 
    xlab("") + ylab("Mean coverage") + theme_bw() + 
    theme(axis.text.x=element_text(angle=60, hjust=0.5, vjust=0.5)) +
    guides(fill=FALSE)
print(qp)
# ggsave(filename = "coverage.pdf", plot = qp, width = 5, height = 5, useDingbats = F)

# Plot coverage for all samples
subtypes = c("pRCC", "chRCC", "RO", "tRCC", "sarcomatoid", "unclassified")
names(subtypes) = c("papillary", "chromophobe", "oncocytoma", "translocation", "sarcomatoid", "unclassified")
col = c(chRCC="#FF3030", RO="#00CC00", pRCC="#800080", tRCC = "#4682B4", sarcomatoid = "#FF8C00", unclassified = "#D3D3D3")

st2 = subtypes[cov_df$Subtype]
col_vec = col[st2]
cov_df2 = data.frame(cov_df, color = col_vec, subtype2 = st2)

get_lim <- function() {
    c(min(cov_df2$cov), max(cov_df2$cov))
}

get_x <- function(x1=0.85, x2=1.15) {
    byy = (x2-x1)/nrow(cov_df2)
    sample(seq(x1, x2, by = byy), nrow(cov_df2))
}

# pdf(file = "coverage2.pdf", width = 5, height = 5)
par(mar = c(4,4,1,1))
boxplot(cov ~ all, vertical = TRUE, outline = FALSE, ylab = "Mean Coverage", ylim = get_lim(), data = cov_df2, xlim = c(0.4, 1.8), xaxt = "n")
points(x = get_x(), y = cov_df2$cov, col = col_vec, pch = 19, cex = 0.7, xaxt = "n")
legend("topright", title = "Subtype", names(col), col = col, pch = 19, cex = 0.6)
mtext(side = 1, text = "Samples", at = 1, line = 1)
# dev.off()


## ======================================================================
## Pathways
## ======================================================================
# This section only produce CSV file
# /gnet/is3/research/data/dnaseq/analysis/thongn/projects/nccRCC_391/pathways_full: calculate how frequent a gene is CNV_gain. CNV_loss, mutated, expression_gain, expression_loss
# /gnet/is3/research/data/dnaseq/analysis/thongn/projects/nccRCC_391/pathways: simular to the pathways_ful but combine all change event (CNV_loss, CNV_gain...) into one (i.e. altered or not)
